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The nonequilibrium reweighting technique, which was recently developed by the present authors, 
is used for the study of the nonequilibrium steady states. The renewed formulation of the nonequlib- 
rium reweighting enables us to use the very efficient multi-spin coding. We apply the nonequilibrium 
reweighting to the driven diffusive lattice gas model. Combining with the dynamical finite-size scal- 
ing theory, we estimate the critical temperature T c and the dynamical exponent z. We also argue 
that this technique has an interesting feature that enables explicit calculation of derivatives of 
thermodynamic quantities without resorting to numerical differences. 
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Most phenomena occurring in nature are in nonequi- 
librium states. Nonequilibrium systems, such as epi- 
demic pj , vehicular traffic , biological network , and 
catalysis 0, have captured a lot of attention. Monte 
Carlo simulation has become a standard tool in scien- 
tific computing, and advanced simulation methods, such 
as cluster algorithms 0, |(| and generalized ensemble 
methods @, H 0, HU, have been developed. However, 
many advanced Monte Carlo methods are not applica- 
ble to nonequilibrium systems. Efficient Monte Carlo 
algorithms for nonequilibrium simulation arc highly de- 
manded. 

Quite recently, the present authors 11] have developed 
a reweighting method for nonequilibrium systems based 
on the Sequential Importance Sampling (SIS) 0, Il3[ . 
With nonequilibrium reweighting, only simulation at a 
single temperature is required to obtain information for 
a range of temperatures. The nonequilibrium reweight- 
ing method differs conceptually from conventional Monte 
Carlo methods. In many Monte Carlo methods, a se- 
quence of micro-states is sampled using the Boltzmann 
distribution. One can interpret this as sampling over a 
"path" generated by the associated Monte Carlo updates. 
Thermodynamic quantities are then averaged over this 
path. In nonequilibrium reweighting, many paths are 
first sampled with a trial distribution that is not neces- 
sarily equal to the Boltzmann distribution. Then ther- 
modynamic quantities are calculated based on the rel- 
ative probability between the trial distribution and the 
Boltzmann distribution. The relative probability is called 
"weights" in literature, which we shall use hereafter. The 
advantage of this is that one could sample many paths 
at one temperature and then calculate required thermo- 
dynamic quantities for a range of temperatures. 

Moreover, Saracco and Albano 0, 0] have proposed 
an effective analysis of nonequilibrium phase transitions, 
in the study of the driven diffusive lattice gas model ^(| , 
using a dynamical finite-size scaling theory. The behav- 
ior of nonequilibrium phase transitions can be extracted 
from short time dynamics [ItI ll8L ITfjl ] . If we combine the 
advantages of dynamical finite-size scaling and nonequi- 
librium reweighting, we can achieve an effective way of 
simulation for nonequilibrium systems. 



In this paper, we apply the nonequilibrium reweight- 
ing metho d II IB to the study of the nonequilibrium steady 
states |2(j, |2l| . We illustrate our method on the driven 
diffusive lattice gas model 16J. We reformulate the 
nonequilibrium reweightin g m ethod and implement very 
efficient multi-spin coding [22t |2j| . We also make modifi- 
cations to the dynamical finite-size scaling relation, which 
was originally proposed by Saracco and Albano [lj, [l5| , 
so that the advantages of reweighting and dynamical 
finite-size scaling can be combined. 

Let us start with explaining the driven diffusive lat- 
tice gas model proposed by Katz, Lebowitz and Spohn 
(KLS) [la]- This system is one of the most well 
known nonequilibrium models exhibiting a nonequilib- 
rium steady state. It was first proposed as a model for 
super-ionic conductors, and attained its popularity due 
to its complex collective behavior. It is constructed as a 
L x x L y square lattice with half-filled lattice sites having 
periodic boundary conditions. Its Hamiltonian is given 
by 

H = -4 nijn V j>, (1) 

where the summation is over nearest lattice sites. The 
variable n,j = 1 when the site is filled and = 
otherwise. Attempts for each particle to jump to an 
empty nearest neighbor site are given by the Metropo- 
lis rate jH], 

T .eWW) = mm [1, exp(-/3(AW - eE))] , (2) 

where a and a' are the system configurations before and 
after the jump, AH. represents the change in energy due 
to the jump, E is a constant driving force, e = — 1, or 1 
depending on whether the jump is against, orthogonal or 
along the direction of the drive, and (3 — 1/T is the in- 
verse temperature of the thermal bath. The L y direction 
is taken as the direction of the drive. The KLS model 
exhibits an order-disorder second order phase transition. 
In the ordered phase, strips of high- and low-density do- 
mains are formed along the direction of the drive. In the 
final steady state, the particles are condensed into a sin- 
gle strip parallel to the direction of the drive [2!| • Hence 
the order parameter can be defined as the density profile 
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along the direction of the drive |bj] , and moments of the 
order parameters are given by 




where n,j = or 1 as defined in Eq. (JIJ, and k = 1,2,4 
represents the first, second and fourth moments of the 
order parameter, respectively. 

We briefly review the nonequilibrium reweighting 
based on SIS, and show the implementation on the KLS 
model. Define a path xt, a sequence of points in phase 
space Oi which were visited during the course of simula- 
tion, as 

X t = (CTl,(72,---<7 t ). (4) 

This path can be sampled using the Monte Carlo method 
at an inverse temperature (3 and a constant drive E. The 
objective is to calculate the appropriate weights for com- 
puting the thermal average of a quantity Q at another 
inverse temperature (3 1 and another drive £", 

n n 

i=i 3=1 

where the sum is over all sampled paths indexed by j 
and w 3 t are the weights. The number of paths is denoted 
by n. To calculate the weights, the following steps are 
implemented, 

1. Suppose x\ — (a{, ■ ■ ■ , a 3 t ) up to time t is sampled 
from a simulation at the inverse temperature (3 and 
drive E. 

2. To go from t, choose a pair of neighboring lattice 
sites at random. If one of the two sites is empty, 
move the particle to the empty site with the rate, 
Ti3,e{&' 3 \&l ), which is the Kawasaki spin exchange 
process. a n denotes the system configuration after 
the move. 

3. An incremental weight Sw J has different values ac- 
cording to two possible outcomes, 

(a) If the move is accepted; of, 1 = a' J and Sw^ = 
T & ,M° n \vl)lT,A° n \vl)- 

(b) If the move is rejected; cr;? +1 = crj and 8w' J — 

ll-Tp,, E ,(a> j \c4)}/ll-Tp, E (<r' j \<4)]- 

The weights at t + 1 is given by this incremental 
weight through the relation w 3 t+1 — 5w J x vxj. with 
w{ = 1. 

For each path j e {I, • ■ ■ ,n}, these steps are repeated 
until some predetermined maximum Monte Carlo time is 
reached. 



We make a comment on the technical detail of calcu- 
lating the weights. For case of infinite drive (E = oo), 
possible values of incremental weights 5wi are, 

Sw = 1, 

Swt = exp(-I2(/3'-/3)), 

Sw 2 = exp(-809' - (3)), 

Sw 3 = exp(-4(/3' - /?)) (6) 

Sw 4 = (1 - exp(-12/3'))/(l - exp(-12/3))), 

Sw 5 = (l-exp(-8/?'))/(l-exp(-8/?))), 

Sw 6 = (l-exp(-4/3'))/(l-exp(-4/3))). 

The weights can then be written as a product of incre- 
mental weights, 

v4 = {5w 1 ) h ^ t \5w 2 ) h3 ^ ■ ■ ■ (6w 6 ) h ° {t) . (7) 

where h\ (t) ■ ■ ■ h J 6 (t) are the number of hits on the in- 
cremental weights Swi ■ ■ ■ Swq during the course of sim- 
ulation from time I to t. Note that Swq is irrelevant 
in Eq. Q). Generalization of this counting method to 
the case of finite E is trivial. Since the calculation of 
weights has been reduced to accumulating a histogram, 
the multi-spin coding technique [2^ can be implemented 
not only for the spin update process but also for the cal- 
culation of histogram of incremental weights. For system 
configuration updates, we follow the multi-spin coding 
technique similar to the case of the Kawasaki spin ex- 
change model [2i| . Once the histogram • • • h^(t) is 
obtained, using Eq. Q) allows us to reweight to a large 
number of temperatures (drives) with negligible extra 
computational efforts. A large increase of efficiency has 
been obtained by a new formulation of the nonequilib- 
rium reweighting. The details of the multi-spin coding for 
the nonequilibrium reweighting will be given elsewhere. 

For the dynamical finite-size scaling, we use the follow- 
ing equation, 

p k = b~* P < k \b~ z T,b^ e,b- l L y ,b~^ L x ,b Xa Po ), (8) 

where k is the fcth moment of the order parameter, p*( fc ) 
is the scaling function for the fcth moment, 6 is the spa- 
tial rescaling factor, e = (T — T c )/T c , (3 is the critical 
exponent for the order parameter (it should not be con- 
fused with the inverse temperature), z/|| and v± are the 
critical exponents for the correlation length parallel and 
orthogonal to the drive, respectively, z is the dynam- 
ical exponent and r is Monte Carlo steps per site. In 
addition to the original scaling form of Saracco and Al- 
bano [Til mi l , our scaling form has a term b x ° pa to reflect 
the initial system configuration 0, ^| . xq is an inde- 
pendent exponent and po <§; I is the order parameter of 
the initial configuration. Simulations have to be started 
with a chosen value of po for all samples. We prepare 
our initial configuration with po = by inserting L y /2 
particles for each column in the lattice and then shuffling 
each column independently. Letting b — r 1 / 2 , we have 

P k =T~^p< k \T^e,T-iL yi T~^L x ,T^ Po ). (9) 
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FIG. 1: Plot of order parameter with infinite drive for 64 x 32 
lattice with actual simulation performed at T = 3.160 shown 
with a bold line. From top to bottom values of T are 3.150, 
3.155, 3.160, 3.165, 3.170. Averages were taken over 4.096 x 
10 6 samples. 



In the limit of L x — > oo at the critical temperature (e = 0) 
with po — 0, the ratio-of- moments reduces to a scaling 
function with a single argument, 



(P 2 ) 2 



= g(r- 1/z L y ) with p Q = 0, e = 0, L x -> oo. (10) 



By plotting the ratio-of-moments versus tL~ z at T c and 
po = 0, neglecting corrections to scaling, the curves 
for different system sizes L y will collapse into a single 
curve. A measure of goodness-of-fit can be defined for 
the "curve-collapse" as 
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\ql v1 { x ) -9L y2 (x)\ dx, (11) 



where g Lyl (x) = g(rL yl z ) and g Ly2 (x) = g{rL y z ). Our 
task is to choose T c and z which minimize r\. In using the 
relation Eq. 11U|) . we should check that the system size 
L x orthogonal to the drive is large enough. At this point, 
we should mention that Leung |2y| has studied the KLS 
model using finite size scaling at noncquilibrium steady 
states. While we focus on dynamical behaviors, his finite 
size scaling method was developed for analysis at stead 
states. 

We now show the results of the Monte Carlo simulation 
for the KLS model. We first illustrate the reweighting 
for the order parameter, and then show how reweight- 
ing can be combined with dynamical finite-size scaling 
(Eq. tHH ) to calculate the critical temperature and dy- 
namical exponent. Figurc^shows how data over a range 
of temperatures can be extracted from simulations at a 
single temperature. The temporal evolution of the order 
parameter p for the infinite drive (E = oo) was investi- 
gated for 64 x 32 lattice. Simulations were performed at 



FIG. 2: Plot of order parameter with finite drive for 32 x 
32 lattice with actual simulation performed at (T, E) = 
(2.765,0.515) and (T, E) = (2.780,0.500). Reweighted data 
are combined using weighted mean. ^From top to bottom 
values of T and E are (T, E) = (2.760,0.520), (2.765,0.515), 
(2.770,0.510), (2.775,0.505), (2.780,0.500), (2.785,0.495). Av- 
erages were taken over 2.048 x 10 6 samples for each simulation. 
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FIG. 3: Plots showing the goodness-of-fit g and correspond- 
ing values of dynamical exponent z for various temperatures. 
Data are generated by fitting ratio-of-moments for L y = 64 
and L y — 128 between the range (p 4 )/(p 2 ) 2 = 1.2 and 
1.405. T c is estimated from the temperature with the best 
fit (T c = 3.175 ± 0.002). Error bars were generated by fitting 
several ranges of ratio-of-moments. 



T = 3.160, and data were reweighted to nearby temper- 
atures, T = 3.150,3.155,3.165,3.170 (from top to bot- 
tom). Averages were taken over 4.096 x 10 6 samples. We 
made independent calculations directly at T — 3.150, for 
example, to check the effectiveness of the reweighting. 
The deviation of the data between the reweighted ones 
from T = 3.160 and the direct ones at T = 3.150 are 
found to be the same within statistical errors. 

We also made simulations for the finite drive (E w 0.5). 
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FIG. 4: Scaling plot of (p 4 )/(p 2 ) 2 versus tL~ z for 2 = 2.09, 
L y = 64 (solid line) and L y = 128 (dotted line) at T = 3.175. 
Initial system configurations were prepared with po = 0. 

We illustrate the reweighting over both E and T. We 
performed two simulations at (T,E) — (2.765,0.515) 
and (2.780,0.500) for 32 x 32 lattice. The reweight- 
ing of the order parameter is made by using p = 
(E l=i Pk I Afc ) / (E L i 1 / A fc ) » where Pi,2 and Ai, 2 are the 
order parameter and error estimates from the first and 
second simulations, respectively. Figure|21shows the tem- 
poral evolution of the order parameter for several temper- 
atures and drives. Data was reweighted to several values 
at (T,E) = (2.760,0.520), (2.770,0.510), (2.775,0.505), 
(2.785,0.495). Averages were taken over 2.048 x 10 6 
samples for each simulation. Generally, we found that 
reweighting is effective when the distributions P ( 9 i _e(xj) 
and Pfji t E'(0{") have sufficient overlaps. Error bars and 
fluctuations of weights [l3l ] can also be used as quantita- 
tive measures on the effective range of reweighting. 

To determine T c , we use the dynamical finite-size scal- 
ing of the ratio-of- moments (Eq. I|10|) ). Here we concen- 
trate on the infinite drive (E = oo). We simulated 64 x 64 
and 64 x 128 lattices, and calculated the ratio of the mo- 
ments, (p A ) I (p 2 ) 2 ■ Before going into the discussion of 
the fitting, we make a comment on the system size L x 
whether we can consider as L x — ► oo. We performed sim- 
ulations for both L x — 64 and L x = 128, and confirmed 
that the ratio of moments for 64 x L y and 128 x L y co- 
incided with each other to within statistical fluctuations. 
Thus, we may regard that L x = 64 is large enough. Since 
v\\ > v±_ for the KLS model, the correlation length or- 
thogonal to the drive, £±, develops slowly; hence L x = 64 
is large enough to use the scaling relation Eq. (|10fl . Now 
we show the fitting procedure. Fitting was performed 
for several temperatures near T c , which were reweighted 
from the data obtained at a single temperature, and for 
each temperature we adjusted the value of z such that the 
goodness-of-fit 77, Eq. Qllfl. becomes minimum. Figure |3| 
shows 77 for several temperatures and the values of z used 
to calculate 77. The best fit occurs at T c = 3.175 ±0.002; 



the error bar on T c is estimated by including all neighbor- 
ing temperatures where the mean values of 77 are within 
two standard deviations of 77 at T — 3.175. The value of 
z within T = 3.175 ± 0.002 is z = 2.09 ± 0.01, and we 
use this value as our estimate of the dynamical exponent. 
Figure^shows the scaling plot of (p 4 )/{p 2 ) 2 as a function 
of tL~ z for 64 x 64 (solid line) and 64 x 128 (dotted line) 
lattice sizes at T = 3.175 and z = 2.09. The curves are 
almost indistinguishable at this scale although some cor- 
rections to scaling can be observed below tL~ 2 = 0.02. 
To study the corrections to scaling, the goodness-of-fit 
for ratio-of-moments for smaller sizes, that is, 64 x 32 
and 64 x 64 lattices, was also calculated using a simi- 
lar procedure. The best fit occurs at T = 3.155 ± 0.005 
with z — 2.23 ± 0.03. The estimate for T c increases with 
the system size, whereas that for z decreases. Our es- 
timates of T c and z are compatible with the recent es- 
timates for infinite lattice, T c = 3.1980 ± 0.0002 [27J, 
T c = 3.200 ± 0.010 15], z = 2.016 ± 0.040 A more 
systematic analysis of the corrections to scaling to get a 
precise estimate of T c and several critical exponents for 
infinite lattice will be left to a separate publication. Be- 
fore closing we show the actual procedure of the reweight- 
ing for each system size. For 64 x 32 lattice, 4.096 x 10 6 
samples were used for the simulation at T = 3.16. For 
64 x 64 lattice, 8.19 x 10 5 samples were used for each 
simulation at T — 3.174 and 3.180. Results were then 
reweighted to other temperatures and combined using 
weighted mean, f = (ELi^/AD/(ELiV^)- Here 
7*1,2 and A12 are the ratio-of-moments and error esti- 
mates from the first and second simulations, respectively. 
For 64 x 128 lattice size, 1.64 x 10 5 samples were used 
for each simulation at T = 3.174,3.177 and 3.180, and 
reweighted results were combined using the same proce- 
dure. 

To summarize, we have studied the use of nonequi- 
librium reweighting based on SIS for the noncquilibrium 
steady states. We have reformulated the nonequilibrium 
reweighting method, which is convenient for the multi- 
spin coding. As a result, a large increase of efficiency 
has been achieved for the performance of simulations. 
We have applied the nonequilibrium reweighting to the 
driven diffusive lattice gas model (the KLS model) . Com- 
bining with dynamical finite-size scaling theory, we have 
estimated T c and the dynamical exponent z. 

Finally, we make a remark on possible applications. 
The nonequilibrium reweighting method is very general 
and has some very interesting properties. For exam- 
ple, the fluctuation-dissipation theorem does not hold 
for nonequilibrium systems and derivatives of thermo- 
dynamic quantities had been estimated using finite dif- 
ferences |28|. With reweighting, derivatives can be cal- 
culated directly by differentiating the weights explicitly, 
that is, 

Y^n ^j\^ w t v-m dw J t 

d P Ej=i w i Ej=i m 
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Here, dw J t /df3' can be obtained by differentiating Eq. Q 
with respect to /?'. We believe that the noncquilibrium 
reweighting method would have several directions for ap- 
plications. 
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